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The phase space for Hamiltonians of two degrees of freedom is usually divided into 
stochastic and integrable components. Even when well into the stochastic regime, 
integrable orbits may surround small stable regions or islands. The effect of these islands 
on the correlation function for the stochastic trajectories is examined. Depending on 
the value of the parameter describing the rotation number for the elliptic fixed point 
at the center of the island, the long-time correlation function may decay as t^^ or 
exponentially, but more commonly it decays much more slowly (roughly as i~^). As 
a consequence these small islands may have a profound effect on the properties of the 
stochastic orbits. In particular, there is evidence that the evolution of a distribution of 
particles is no longer governed by a diffusion equation. 
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1. Introduction 

Many important problems in physics are described by Hamiltonians of two degrees of 
freedom. Examples are the motion of a charged particle in electrostatic waves, the motion 
of a charged particle in various magnetic confinement devices, the acceleration of a particle 
bouncing between a fixed and an oscillating wall, the wandering of magnetic field lines, etc. 
In such systems, there may be a parameter k which governs the strength of the coupling 
between the two degrees of freedom. When fc is zero (no coupling), the system is integrable. 
As k is increased, some of the integrable trajectories disappear, and the motion becomes 
stochastic. Eventually, when the k is large, nearly all of phase space is occupied by stochastic 
trajectories. These properties are nicely illustrated by the standard mapping [1] . Although 
this mapping is an idealization, it accurately portrays many of the properties of real systems. 

When k is small, various perturbation theories are available to describe the trajectories 
approximately. On the other hand, in the highly stochastic regime (k large) certain statis- 
tical quantities may be analytically determined by assuming that the motion is ergodic and 
that the correlation time is short. In the case of the standard mapping, this allows simple 
determinations of quantities such as the diffusion coefficient [1, 2] and the KS-entropy [1]. 

The purpose of this paper is to examine more critically the assumption that the corre- 
lation time is short in the stochastic regime. Our interest in this problem was triggered by 
studies of the diffusion coefficient for the standard mapping, 

- n-i = -fcsin6'i_i, 0t-0t-i=rt. (1) 

Far above the stochasticity threshold fc ^ 1, the diffusion coefficient defined as 

P=lim«!l^, 
t—^oo 2t 

where the average is over some appropriate ensemble, is given by assuming that ^ is a random 
variable in the equation for r. This is equivalent to assuming that the correlation function is 
proportional to a delta function and gives the "quasi- linear" rc;sult V/Vqi = jk"^. Including 
the correlations out to short times (about t = 4) gives corrections to the diffusion coefficient 
reported by Rechester and White [2] which can enhance the diffusion coefficient by as much 
as a factor of about two. However, a numerical determination [3] of the diffusion at fc = 6.6 
gave V/Vqi ~ 80. At this value of k there is an island ("accelerator mode") present in the 
stochastic sea. Although the orbits used to compute T) were all in the stochastic region, they 
were able to wander close to the island and stay close to it for a long time. This introduces 
long-time correlations into the motion and accounts for the anomalously large value of V 
observed. 

The appearance of islands in the stochastic regime is not at all unusual. This may 
be seen from Sinai's estimate [1] for the munber viT) of periodic orbits of period < T: 
v{T) ~ exp(/ir) for T oo where h is the KS-cntropy. Now the majority of these periodic 
orbits are unstable. However as k is increased, h generally increases {h ^ log(ifc) [1] for 
the standard mapping) and hence new periodic orbits must appear. Generally a tangent 
bifurcation is responsible for the appearance of the new periodic orbits. At the tangent 
bifurcation, a pair of stable and unstable periodic orbits is born. The stable orbit gives rise 
to other longer periodic orbits as k is increased and eventually becomes unstable. Thus as 
we increase k, many small islands appear via tangent bifurcations, survive for some interval 
of k, and finally disappear (usually through period doubling). If we pick a particular value 
of fc, it is not clear that there will necessarily be any islands present. However, we may 
speculate that at some arbitrarily close value of fc, there will be some islands. 
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It may be objected that the large effect seen in the standard mapping arises because 
the islands are accelerator modes [1] and that such islands are a rather special feature of 
the standard mapping. While accelerator modes are the only islands which will contribute 
significantly to the force (i.e., acceleration) correlation function and hence to the diffusion 
coefficient, any islands will contribute to the correlation of some functions on phase space. 
The results of this study will be applicable to systems like the Fermi map which have no 
accelerator modes. This study also contributes to the understanding of the more general 
problem of motion in a divided phase space. 

In this paper, we wish to examine more closely the effect these islands have on a 
stochastic trajectory. As far as determining the effect on the correlation function, this 
involves determining how "sticky" the island is. Given that the stochastic trajectory comes 
within a certain distance of the boundary of the island, how long do we expect it to stay 
close to the island? This approach is inspired by work of Channon and Lebowitz [4] on the 
correlations of a trajectory in the stochastic band trapped between two KAM surfaces in 
the Henon map. Similar work has been carried out on the whisker map by Chirikov and 
Shepelyansky [5]. This work is being extended by B. V. Chirikov and F. Vivaldi. 

Since we concentrate only on the behavior close to the island, this approach may be 
characterized as a local one. This should be compared with Fourier transform methods 
[2], which are global and are not well suited to the description of localized phenomena. 
For instance, Meiss et al. [6] attempted to use such methods to compute the long-time 
correlations for the standard mapping, and they found poor agreement with numerical 
experiments whenever islands were present. 

The paper is organized as follows. In section 2, we derive a canonical mapping which 
describes the behavior near an island. Next (section 3), we define the trapping statistics 
which describe how sticky the island is. The results for the trapping statistics are given 
in sections 4 and 5. In section 6, we show how to apply these results to obtaining the 
correlation function. The results are discussed in section 7. 



2. Derivation of mapping. 

Far into the stochastic regime for a general mapping, the islands which appear via 
tangent bifurcations are very small and exist only for a small interval in parameter space. 
This allows us to approximate them by a Taylor expansion in both phase and parameter 
space about the tangent bifurcation point retaining only the leading terms. This was carried 
out in ref. 3 where the resulting mapping was reduced to a canonical form 

Q : Ut- yt-i = 2(a:t_i - K) = g{xt-i;K), Xt - xt-i = yt- (2) 

Here K is proportional to /c — fctang (fctang is the parameter value where the tangent bifurca- 
tion take place) and x and y are related to the original phase space coordinates by a smooth 
transformation. The mapping Q represents an approximation of the general mapping close 
to the point of tangent bifurcation. For K < I, Q has no periodic orbits. At K = 0, there 
is a tangent bifurcation when an unstable fixed point appears at x = y = 0. (This is not a 
hyperbolic point since its stability is determined by the quadratic terms in the mapping.) 
For < K < 1, this fixed point splits into a pair of stable (elliptic) and unstable (hyper- 
bolic) fixed points located at {x,y) = (=FV^, 0), respectively. The elliptic fixed point is 
usually surrounded by integrable trajectories (KAM curves) which define a stable region 
(the island) in which the motion is bounded. An example of island structure is shown in fig. 
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1 for if = 0.1 (the value of K at which extensive numeric;al cialculations have earried out). 
At if = 1 the stable fixed point becomes unstable and gives rise to a period-2 orbit via a 
period-doubling bifurcation. At if = 1.2840 the period-doubling sequence accumulates [7], 
and at this point (or shortly hereafter [8]) the area of the stable regions becomes very small. 
For < if < 1, the mapping Q may be transformed to the Henon quadratic map [9], the 
parameter if being related to Henon's cos a by 



{a is the mean angle of rotation for points close to the stable fixed point). 

Referring again to the islands shown in fig. 1, consider a particle which at t = is 

close to, but outside, the islands. (We often speak of an orbit in terms of the position of a 
particle whose equation of motion is given by Q.) Initially, the particle will stay close to the 
islands; however as we let t — > ±oo, we find {x,y) — > (oo,±oo). It is just such trajectories 
we are interested in, because they correspond to particles in the stochastic region of the 
general mapping approaching the islands, staying there for some time (and contributing to 
long-time correlations), and then escaping back to the main part of the stochastic region. 

What we want to do is to follow such trajectories numerically and see how long they stay 
close to the islands. However, we need some method of fairly sampling these trajectories. 
"Fairly" means that we should sample them in the same way that a stochastic trajectory 
of a general mapping does. Since the stochastic trajectory is crgodic over the connected 
stochastic region of phase space, we must sample trajectories in the same way; i.e., we must 
ensure that the superposition of all the sampled trajectories covers the region outside the 
islands uniformly. 

We achieve this by changing Q so that the phase space is compact. This may be 
accomplished by replacing g{x;K) in Q by the periodic function 



where L = Xmax — iCmin > 0. The resulting map will be called Q* . If Xmin and Xmax are 
chosen to span the region where there are islands in Q, then Q* obviously will contain the 
same islands. Furthermore the motion close to the islands will be the same. By replacing 
g by g* , the whole of phase space becomes periodic in the x and y directions with period 
L. The motion can be treated as though it were on a torus. A particle which starts near 
the island will, as with Q, spend some time close to the island. But when it moves away 
from the island it no longer goes to infinity, but rather it loops around the torus and has 
another chance to approac;li the islands. Since Q* is area-presc^rving, this single orbit will 
ergodically cover the region outside the islands in the desired manner. Examples of such 
orbits are shown in fig. 2 for the same parameters as for fig. 1. (One embarassing feature 
of Q* is that new islands are introduced. They are in fact accelerator modes a particle 
inside one of them loops around the torus in either the x ov y directions. As discussed in 
appendix A, these islands do not effect the statistics for the long trapping times.) 

One useful way of looking at Q* is as a magnification of a small region near a tangent 
bifurcation in the general mapping. The difference is that once the trajectory leaves the 
vicinity of the islands, it is immediately re-injected on the other side of the islands. In the 
general map, the trajectory will spend some long time, which depends on the ratio of the 
size of the islands to the total accessible portion of phase space, in the stochastic sea before 
coming back to the vicinity of the islands. 



cos a = 1 — 2Vk 




otherwise. 



< X < x^ 



max; 



(3) 
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Assuming that the long-time behavior of stochastic orbits is dominated by the region 
close to the islands, there are two advantages to reducing the problem to a study of Q* . 
Firstly, since Q* describes the behavior of most islands far into the stochastic regime, the 
properties of many mappings may be treated by looking at a special mapping Q* which 
depends only on a single parameter K. The second advantage is that the properties of 
orbits close to the islands may be studied much more efficiently because there is no need to 
follow orbits while they spend a long and uninteresting time far from the islands. 



3. The trapping statistics. 

The prescription for numerically determining how sticky is the island system in Q for 
a given K is to pick trajectories outside the islands in Q* and to iterate the mapping many 
times. The mapping is performed on the torus, i.e., x and y are reduced to their base 
intervals [xmin, Xmax) and [—\L,^L) after each iteration. However, we keep track of when 
an orbit moves off one edge and re-appears at the opposite edge. The orbit is then divided 
at those points when the orbit looped around the torus, and the lengths of the resulting 
orbit segments are recorded. The main result of such a calculation are then the trapping 
statistics ft which are proportional to the number of orbit segments which have a length of 
t. 

Suppose that the total length of the orbit is T and Nt is the number of segments of 
length t. If T is so large that we can ignore partial segments at the ends of the orbit (this 
problem is examined below), then we have Y^tNt — T; the total number of segments is 
N = Yli^t- The trapping statistics are defined by ft = N^/T and are therefore normalized 
so that J2^ft = 1- The mean length of the orbits is given by a = (= T/N). The 

probability that a particular segment has length t is pt = aft {— Nt/N). If an arbitrary 
point is chosen in the orbit, then tft is the probability that this point belongs to a segment 
of length t and ft is the probability that it belongs to the beginning, say, of a segment of 
length t. 

Three factors effect the measurement of ft- They are (a) the presence of the spurious 
accelerator modes, (b) the choice of a::niin and a:max, and (c) the total length T of the 
trajectory used to measure ft- The first two items only effect ft for small t (apart from an 
overall normalization). In order to account for the last item, we define ft hy Nt/{T +1 — t) 
(rather than Nt/T). This accounts for the fact that we are less likely to observe orbit 
segments whose length is close to T. All these points are discussed in detail in appendix A. 

The survival probability 

oo 

Pt=Y.Pr (4) 

T=t + 1 

is the probability that an orbit beginning in a segment at t = is still trapped in the same 
segment at time t. Note that Pq = 1 as required. This is the quantity studied in refs. 4 and 
5. The correlation function 

oo oo 
Cr = Y.^t-T)ft=Y,Pt/0^ (5) 

is the probability that a particle is trapped in the same segment at two times r apart. Again, 
we have Co = 1. 
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Thcrc arc two other ways of interpreting Cr- If we start many particles at positions 
uniformly distributed in the stochastic sea of Q* (i.e., in the dark region of fig. 2), then Cr 
gives the fraction of particles remaining in the L x L square after r iterations of Q (rather 
than Q*). Alternatively, consider a drunkard who executes a one-dimensional random walk 
with velocity v = dr/dt = ±1. The direction of each step is chosen randomly, while the 
durations of the steps are chosen to be the lengths of consecutive trapped segments of Q* . 
Then for integer r, Ct is just the usual correlation fimction for such a process, i.e., (wtf t+r)<- 
The behavior of this random-walk process is similar to the behavior of an orbit in the general 
mapping when two accelerator modes with opposite values of the acceleration are present. 
(This is the case with the first-order accelerator modes for the standard mapping.) In section 
6, we will show how Cr is related to the correlation function for the general mapping. 

A diffusion coefficient may be defined by 



This gives the diffusion rate for the drunkard in the random- walk problem above. It is also 
related to the diffusion coefficient for the general mapping (see section 6). 

Since Q* must be iterated many times to provide good statistics for /( for large t, 
extraordinary steps were taken to ensure that the numerical program was fast and reliable. 
The time for one iteration on a Cray-1 is a little less than 75 ns. One way that the code was 
made more reliable was by doing the arithmetic in fixed-point (as opposed to floating-point) 
notation. The numerical mapping is then one-to-one which is the discrete counterpart of 
area-preserving. (Floating-point realizations of mappings are typically many-to-one.) This 
precludes the possibility of an orbit, which starts far from the island, approaching the 
island and becoming permanently trapped near the island. Even though such behavior is 
forbidden for an area-preserving mapping, such behavior may be observing with a floating- 
point realization. Details of the numerical methods are given in appendix B. 



4. The results for small K. 

We begin by considering the cases where K is small or zero. In this case the mapping 
equations are nearly integrable and this enables us to derive approximate analytic expres- 
sions for ft- Figure 3 shows ft for K = —10"'', 0, and 10"''. Three types of behavior are 
seen for t ^ oo: a cutoff distribution ft = for t > fmax ~ 50, a rapid algebraic decay 
ft ^ t~'^, and an exponential decay ft ~ exp(— O.lt). 

The easiest case to begin with is K = 0. The method for deriving ft analytically 
consists of computing the length of the trajectory through some point and then assigning 
some probability that this trajectory will be chosen. The first part of the calculation has 
been carried out by Zisook [10] . We repeat the calculation here to establish the method for 

other cases. 

When K — 0, we are exactly at the tangent bifurcation point. There are no islands in 
this case, but trajectories can still spend arbitrarily long near the fixed point at (x, y) = 
(0, 0). Near this point Xt — Xt-i and yt — yt-i are small. We therefore rcscalc x, y, and t 
with X = e°'X, y = e^Y, t = e~^T where e is small. The mapping Q becomes 



oo 




(6) 



T = l 
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X{T) - X{T - e) 



= e 



^Y{T). 



e 

Choosing a = 3 and /? = 2 and replacing the left hand sides by derivatives, we obtain 

These are Hamilton's equations (with X and Y being conjugate position and momentum 
coordinates) for the Hamiltonian 

Curves of constant H in the {X, Y) plane give the trajectories, examples of which are shown 
in fig. 4(a). We define the trapping time as the time it takes to traverse one of these curves 
from Y = — oo to oo. (This Hamiltonian gives escape to infinity in a finite time. The 
time taken for a particle to escape to infinity in Q is infinite, but very weakly so. The 
particle roaches y from y = in roughly log logy steps for y large.) Since dX/dT = Y = 

^2{H + ^X^), the trapping time may be written as 



T(Xo) =2 [ 



dX 



Xo ^|(X3-X3) 

where Xq, = — (|iJ)^/"^ is the X intercept of the trajectory with Y = 0. Performing the 
integration gives 

T(Xo) = X 4.207 \Xo\-^'^ , for X ^ 0. 

(The numbers here may be written in terms of incomplete elliptic integrals.) 

This completes the computation of the trapping time. We now assign weights to each 
trapping time by requiring that particles spend equal times in equal areas of phase space. 
Let A(X[)) be the area between the trajectory passing through the origin and that passing 
through (Xo,0) in fig. 4(a). Since Y ~ X^/"^, scaling invariance gives A{Xq) = A(\)X'^^'^ . 
(This procedure needs to be carried out separately for positive and negative Xq. However, 
the scaling relations are the same in the two cases.) Parameterizing in terms of the trapping 
time T gives A{T) ~ T~^. The fraction of particles which are trapped for times between 
T and T + dT is proportional to the differential area dA{T) ~ T~^dT. Finally, we divide 
by T to give T~'^dT as the number of orbit segments of lengths in this range. In unsealed 
variables, we have ft ~ t~'^ which is valid for t large. The correlation function has the 
asymptotic form Cr ^ . 

It is interesting to enquire what the asymptotic behavior for ft would be if g in (2) 
were a higher-order polynomial in x. If g{x; 0) = x™, with m > 1, we can repeat the above 
calculation and find that 

. 1 

^(3m+l)/(m-l) ■ 

Since /j has a finite second moment, D always exists. (For to 1, we find an exponential 
decay of ft- This corresponds to the case K > discussed below.) 
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Whon K is small and negative, there are no periodic orbits. The particle; can spcind 
only a bounded time close to the origin. Defining scaled variables as before, together with 
K = —e^, gives differential equations which are derivable from the Hamiltonian 

H=\Y'^- - 2X. (8) 

The trapping time is now 

T(Xo) = r (9) 
Jxo ^lx^ + X-\Xl-Xo 

where Xq is the X intercept of the trajectory with F = 0, i.e., it is the real root of H+ |-^o + 
2Xa = 0. This is plotted as a function of Xq in fig. 5(a). We see it attains a maximum value 
of Tinax = 5.1454 at X^ = —0.5536. In unsealed variables this means that the maximum 
trapping time is tmax = 5.1454 \K\ and that this trapping time is attained by particles 



passing through {x,y) = (— 0.5536 0). The \K\ scaling of the maximum trapping 
time has been derived by Zisook [10]. 

To assign probabilities to the various trapping times, we define A{Xo) as the area 
between the trajectory passing through (0,0) and that passing through {Xo,0). Since Y = 

2^iX3 + X- - Xo, we obtain 

A{Xo) = 4 / ^J\X^ + XdX - 4 / ^iX3 + X - - Xo dX. 

(The two integrals need to be done to together to get a finite answer.) Differentiating gives 

dA{Xo)/dXo = 2{X^ + l)T(Xo). 
The number of orbit segments of length T is then proportional to 



dTiXo) 



dXo 



E 



2(^g + l) 
\T'iXo)\ • 



(10) 



The right hand side is written as a function of Xq. This is converted to a function of T 
by inverting T{Xo). Since this gives a double- valued function, the two branches must be 
summed over as indicated by the summation sign. In unsealed variables we have 

for t = Odii'l"^/^). The function F{T) is plotted in fig. 5(b). For T <C T^ax, we have 
F{T) w 6.205 X lO'^T-^ while for T « T,„ax, F{T) w 3.024(rmax - T)"^/^. The correlation 
function is given by the second integration of /* so that for r « tmax we have Cr ~ (^max — 

r)3/2 for T < t^ax. 

Figure 5(b) should be compared with fig. 3(a). The effect of the singularity in F{T) 
at Tmax is evident, although its position isn't quite right. Furthermore, the decay of ft for 
smaller times is somewhat slower (approximately as t~^) than for F{T). These discrepancies 
arise because we are not far enough into the asymptotic regime since e is not very small 
(0.1). The formula t^ax = 5.1454 \K\ may easily be verified for smaller values of \K\. 
The T-'^ behavior of F{T) has of course the same origin as that of ft for K = 0. However 
the numerical results for K = given in fig. 3(b) show that it is not attained until about 
t ~ 100. When K = — 10~^, tmax is only about 50, and there is no interval in which the t~'^ 
behavior is exhibited. 
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Finally wc turn to the case K > Q. As with ii' < 0, we can approximately derive the 
motion from the Hamiltonian 

H = -IX^ + 2X, (11) 

where we have used the same scaled variables as previously except that K = e'^. The 
trajectories for this Hamiltonian are given in fig. 4(b). There is single island centered at 
the fixed point at {X,Y) = (—1,0) and the island extends all the way to the separatrix 
emanating from the unstable fixed point at {X,Y) = (1,0). This is an idealization because 
for the mapping there is a stochastic band close to the separatrix. However, for small K 
this band is very thin, and the island does have quite a "clean" outer edge. 

We may carry out the analysis given for K < Q with appropriate changes to obtain 
/( analytically. However, we may avoid doing a lot of tedious integrals by concentrating 
only on the long-time behavior of ft- For shorter times, namely for 100 < t <^ K~^/^, we 
expect that ft ^ because the relevant trajectories never get close to the island and the 
Hamiltonian (11) approximately reduces to that for = (7). Since only the region close 
to the unstable fixed point contributes to ft for large t, we need only consider this region. In 
addition, we should distinguish those trajectories which encircle the island and so encounter 
the region near the fixed point twice from those which do not and only encounter this region 
once. For a given distance from the fixed point, the trapping time of orbits in the former 
category will be about twice as long as those in the latter category. 

Near any hyperbolic point there exist coordinates (^, r]) in which the mapping may be 
written as 

where A is the Lyapunov number at the fixed point. The trajectories lie on hyperbolae of 
the form S^r] = ±^q. Since the island encircling trajectories dominate the long-time behavior 
of ft, we define the trapping time t{^o) as the twice the time it take to get from t] = 1 to 
^ = 1. (The factor of two accounts for the two encounters with the fixed point.) We find 

t{^o) = -41og^o/logA, dt{^o)/d^o = -4/(CologA). 

The area A{^o) between the hyperbola, the axes, and the lines ^ = 1 and = 1 is 

^(^o)=^o(l-21og^o). 
The trapping statistics are then given by 

1 dA 

.ft ~ -^-^ ~ exp(-ilogAt). 
For the fixed point at (0. \/K), 

X=l + 2Vk + 2\Is/K + K « 1 + 2K^''^, ft ~ exp{-K^'H). 

Cr behaves in the same way. For K = 10~^, the decay rate should be about 0.1, which is 

indeed what was observed in fig. 3(c). Similar agreement is seen aX K = 10~^ and 10"^. 
However at K = 0.1, the central island has shed a chain of sixth-order islands and the 
foregoing analysis does not apply. This case is discussed in the next section. 
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5. The results for K = 0.1. 

Wc have measured fi for K between and 1.3 at intervals of 0.05, and at most of the 
values oi K a, slow algebraic decay of fi is seen. A representative case is if = 0.1, whose 
trapping statistics are given in fig. 6(a), which illustrates the slow decay for very long times 
t ~ 10''. Also given in fig. 6 are Pt, C,-, and a = — dlog C^/dlogr (thus locally Cr ~ t~"). 
This last plot shows the power at which C^- decays varying between about j and |. 

A glance at fig. 2 shows the origin of this behavior. The central island is surrounded 
by a chain of sixth-order islands. Around each of these islands are several other sets of 
islands. This picture repeats itself at deeper and deeper levels. A particle which manages 
to penetrate into this maze can get stuck in it for a long time. 

For T < lO'', fig. 6(d) gives a fa j. Correspondingly we have Pt ~ where p 
■ 5/4. This is close to the asymptotic (t — > oo) result found in ref. 5 for the whisker 
map, in which (p) « 1.45. This is another indication that the behavior of a Hamiltonian 
with a divided phase space has "universal" properties. However, in our case, a shows some 
strong variations beyond t » IC* where Cr "steps down" (e.g., between 10"* and 3 x 10^). 
This means that the asymptotic form of Ct is very difficult to determine numerically. 

The diffusion coefficient D is given by the summation of C and is approximately 6400 ± 
800. The error is estimated by calculating D separately for subsets of the orbits sampled. 
Unfortunately, because a few very long orbit segments have such a large effect on D, the 
individual observations of D come from a highly skewed distribution and the error may be 
severely underestimated. We will try to get a idea of the error by asking what behavior is 
possible for /( for ia = 10* < i < 2 x 10^ = h {h is the length of the orbits used to compute 
/( in fig. 6). Since no segments were observed in this range, we have 



where M = 1600 is the number of orbits used. This just says that the expected number 
of orbit segments in this range of t is less than 1 (see appendix A). Suppose that /( = 
^(i/^a)"^^"*""^ where a is the exponent at which Cr decays and A is the value of /( at 
t = ta = 10®. For tb '> ta and a > 0, we can evaluate the integral to give MAtatb/{'^ + ct) 
approximately. If we take A = 10~^° (this value was estimated from fig. 6(a)), we find 
a > 2. In the case of the slowest decay, a = 2, the portion of ft between ta and oo would 
increase D by about a factor of 2 over the value given above. If ft takes another step down 
near t = 10*, then A might be smaller and smaller values of a would be possible and the 
maximum error in D would be larger. For instance with A = 0.3 x lO^'^'^, then all values 
of a < are consistent with the numerical observations. Since Cr sums to infinity for all 
a < I, D may well be infinite! 

If D is indeed infinite, we would wish to know how a group of particles spreads with 
time. We again consider the drunkard's walk based on Q* which was introduced in section 
3. The second moment of r is related to the correlation function by 



This is plotted in fig. 7(a), using the data of fig. 6. For t < 10^, Si grows somewhat faster 
than (see fig. 7(b)) and even until t « 10^ St is growing significantly faster than 
linearly. Beyond 10^, the numerical data shows a convergence to a linear rate; but this is 





(12) 



r=l 
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merely because no segments longer than about 6 x 10^ were observed. For t ^ oo, St grows 
as assuming that the exponent a at which d- decays asymptotically is less than 1. If 

the diffusion coefficient is estimated from = \St/t, then grows with t as shown in fig. 
7(c). 

When applying these results to the general mapping, we will need to know the area 
&2 occupied by the stochastic trajectories (the dark area in fig. 2). This was calculated by 
dividing the phase space into 1024 x 1024 little boxes, iterating the mapping many times, 
and counting the number of occupied boxes. It is important iterate the map enough times 
so that (a) the expected occupation number of each box is reasonably large and (b) the 
trajectories have time to wander into all the nooks and crannies. (In practice, the second 
requirement is more stringent.) With Xmin, a^max, and Xj. as given in the caption to fig. 6, 
the area of the stochastic component is found to be about 62 = 0.693. 

6. Application of the results. 

We wish now to determine the contribution of an island to the correlation function 
of an orbit in the stochastic component of phase space of a general two-dimensional area- 
preserving mapping G. (The analysis applies equally well to Hamiltonians with two degrees 
of freedom. The phase space is then the Poincare surface of section and the unit of time is 
the period of the island.) 

For simplicity we begin by considering the case where there is a single small island em- 
bedded in the connected stochastic component. Let the total area occupied by the stochastic 
component be Ax . Suppose a small island centered at xq is immersed in the stochastic sea. 
When iterating Q* , we have been approximating G in a small region around xq. The square 
defined by aJmin and Xmax in Q* is transformed into a small box (parallelogram) of area Bq. 
The ratio of the areas in the two spaces is Bo/L"^ = 7. Suppose the area of the stochastic 
component of G which lies inside this box is Bi. (In the notation of appendix A, Bi = 761.) 
Recall that Q* also contains spurious islands (accelerator modes) which have no counter- 
part in G. To account for these islands we define as in (Al). From this we can derive 
^/j* = 1/a* and pi = a* (parallelling the definitions made in section 3). 

It is useful to begin by forming an idea of what a stochastic trajectory will look like. 
Orbits will consist of alternating trapped and free segments. The trapped segments are those 
which are restricted to Bi , while the free segments are those excluded from Bi . (Here and in 
the following we use Bi, etc., to refer to a particular subset of phase space as well as the area 
of this subset.) The basic assumption is that each visit to the island is uncorrelatcd with the 
previous one. So, on first entering the area Bi, we assume that a segment of length t will 
be chosen randomly with probability p* . A simple model for the motion in the stochastic 
region Ai — Bi, which excludes the region near the island, is as follows. The first point after 
a trapped segment is randomly (and with a uniform distribution) situated in the Ai — Bi. 
This is in accord with the picture that once an orbit leaves Bi. it rushes away from the 
vicinity of Bi extremely quickly. If this point is the pre- image under G of _Bi , then the next 
point is the first point of another trapped segment. Otherwise another point is chosen at 
random in Ai — Bi and the procedure is repeated. The mean length of trajectories trapped 
in Bi is a* . The area of the points which are initial points of trapped segments is therefore 
Bi/a*. These are the points whose pre- images lie outside Bi. Thus the probability that a 
point in Ai — _Bi is a pre-image of Si is e = (Bi/a*)/{Ai — Bi). The probability that a 
particular free segment has length t is then e(l — e)*^^. These probabilities ensure that ratio 
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of time that the trajectory spends in Bi and in Ai — Bi is in the ratio of the area of these 
regions. (The assumptions made to obtain the distribution of lengths of the free segments 
is probably overly restrictive. However such a "memory-less" model is probably accurate 
for the long times we are interested in.) This model is discussed in more detail in appendix 
C where it is used as a basis for constructing a Markov-chain approximation of the motion. 

Consider the correlation function 

C(r) = (Mx(t))Mx(t + T)))t, (13) 

where h is some smooth function of the position in phase space x (in particular we require 
that it is a constant throughout Bi). We shall assume that the mean value of /i(x(t)) is zero 
for the stochastic orbits. We identify those terms in (13) for which x{t) and x(t + t) belong 
to the same trapped segment as Cis(r) the contribution to C(t) due to the island. (Thus for 
such terms we have x{t + t') G Bi for all r' such that < r' < r.) Except for r = 0, the 
other terms arc smallcir by a factor of about Bi/Ai, which is typically very small. This is 
so because h has a zero mean and because of the rapid mixing of orbits in Ai — Bi. This 
question is examined in appendix C where it is also shown that the additional terms do not 
contribute to the diffusion coefficient. Cis(T) may be written as 

Ci,(r) = /z2(xo)-^^(i-r)/;. 

The first factor arises because both endpoints are in Bi and h{x) /i(xo) for such points. 

The second factor is the probability that x{t) lies in Bi, and the sum is the probability 
that x{t + t) belongs to the same trapped segment as x(i). This sum is just the correlation 
function defined in terms of /(* instead of ft- For large t, we can substitute for /j* using 
(A2), and the sum becomes {bi/bijCr- Thus the contribution to the correlation function 
due to the island is 

Cis(T) =/i'(xo)^62a. 

Equation (A3) shows that this result is independent (for large r) of the choice of a;max and 
Xrain (which is as it should be). 

If h is the rate of change of one of the components of x, e.g., /i(x) = dx/dt, then the 
island enhances the a;-space diffusion coefficient by 

Di, = /i2(xo)^62A 

where D (assuming that it exists) is given in (6). Whether or not D exists, the mean squared 
X position of a group of particles initially concentrated in a small region is enhanced by 

where Sl is given by (12). 

If there are more than one island, then the contributions should be added together in 
Cr and V. If there is a chain of A^th order islands at xq, xi, . . . , xjv = xq, their contribution 
to the correlation function C{Nt + j) is 
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whcrc 7 = So/ for one of the islands and < j < N. The contribution to V (assuming 
again that it exists) is 



and similarly for Sis {t) 

Let us apply these results to the first-order accelerator modes for the standard mapping 
(1). These are first-order islands which appear at fc = 27rn (n an integer) [1]. The accelera- 
tion in such a mode is r* — rt-i = ±2im. These modes appear in pairs, one with each sign 
of the acceleration. We will take the area of the stochastic component Ai to be equal to the 
entire area of phase space (27r)^. The relation between the parameters is found by matching 
the residues at the stable fixed point. This gives k'^ = (27rn)^ -|- IQK. When transforming 
to the quadratic map Q, lengths are magnified by a factor ^im [3] and so 7 = (2/7rn)^. 
If we wish to estimate the diffusion coefficient, we must define h to be the acceleration; 
thus /i(xo) = ±27rn. For K = 0.1, we take 62 = 0.693 and the proportionality constant 
connecting the scripted and unscripted quantities above is 2/i^(xo)(7/j4i)62 ~ 0.56. (The 
factor of 2 accounts for the presence of the two islands.) 

Using Vq\ = {-Kuf and taking 6400 as a lower bound for D, we find that the contribution 
to the diffusion coefficient is increased over its quasi-linear value by a factor of at least 
360/n^. Thus for n = 1 or fc « 6.41, the islands completely dominate the diffusion. The 
first-order accelerator modes continue to have such a large effect at least until k ~ 100. If 
D is in fact infinite, even arbitrarily small accelerator modes will eventually dominate the 
motion and fig. 7 can be used to estimate the time at which the accelerator modes become 
important. 

7. Discussion. 

We have looked at the effect of a small island on the correlation function for the stochas- 
tic trajectories of Hamiltonians of two degrees of freedom. When the parameter K is small, 
the problem may be treated analytically and we find that the contribution to the correlation 
function is zero for r > 5.1454 when K < Q, decays as r^^ when if = 0, and 

decays as exp(— ii'~^/^r) for K > Q. 

The more interesting case is when K is not small and the island is surrounded by other 
islands. In the case we considered in detail K = 0.1, the decay of the correlation function 
is algebraic and very slow (roughly as t~^) out to times on the order of r ~ lO''. Even 
when the islands are small, this can still lead to an enormous increases of quantities such 
as the diffusion coefficient. Although it has not been definitely established, there are strong 
indications that the diffusion coefficient may be infinite, indicating that the distribution of 
particles docs not obey a diffusion equation and that the particles spread more rapidly than 
diffusively. Such behavior is fairly typical, having been observed at several other values of 



It would be interesting to know how the distribution of particles docs evolve in time. If 
we consider a system like the standard mapping at a parameter value for which an accelerator 
mode exists, then for large t this distribution could, in principle, be found from ft, but its 
determination is beyond the scope of this work. For now, we observe that the distribution 
will be far from Gaussian for times at least until t = 10®. It will contain a very small 
but extremely long tail that contributes significantly to the variance. This has important 
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consequences for minierical experiments. Imagine trying to measure Sc for t = 10^ by 
directly measuring rj — ro for N particles. If N is merely some "reasonably" large number 
like 1000 then St will most likely be greatly underestimated. (To counteract this there is a 
tiny probability that St will be fabulously overestimated.) We know that we should sample 
some orbits with rt — ro ~ 10'' to be able to estimate St accurately. But since St ~ 10^, we 
need to sample at least (10*^)^/10^ = 10^° orbits before the effect of such a long segment 
is correctly diluted. Obviously such a calculation is totally out of the question. A much 
better approach is to measure the correlation function and to derive St using (12). On the 
other hand, the requisite number of orbits probably are sampled in real experiments. For 
instance in plasma physics applications the total number of orbits is typically 10^^. 

There are several questions still to be answered. What is the long-time behavior of the 
correlation function? Figure 6 shows that it has not attained any well-defined asymptotic 
limit by r = 10^. What determines the asymptotic behavior of the correlation function? The 
simplest picture we can form for treating this problem would go something like this: There 
is some outer KAM curve marking the boundary of the main island centered at {—Vk, 0). 
According to Greene [11], this curve is approximated from the outside by a sequence of 
islands whose winding numbers are the rational approximants to the irrational winding 
number of the KAM curve. Thus the long time behavior of ft may be found by considering 
how an orbit wanders through these islands to approach the KAM curve. Greene [12] found 
an algebraic decay of the correlation function based on a simple model of such a process. 
Equivalently this behavior may be obtained from a diffusion equation in which the diffusion 
coefficient approaches zero sufficiently fast as the KAM curve is approached. This picture 
was the one proposed in ref. 5. 

Such pictures are however probably incomplete. There is no reason to suppose that 
the KAM curve around the central island determines the asymptotic behavior of ft- For 
K = 0.1, it could equally well be the last KAM curve around the sixth-order islands which 
surround the central island, or the KAM curve around one of the chains of islands around 
the sixth-order islands, or all of them together! For instance, the longest orbit segment seen 
for K = 0.1, whose length was about 6 x 10^, spends most of its time around such a chain of 
islands whose order is 6 x 23 = 138. (One of the members of this chain is visible in fig. 2(b) 
a,t x + \fK K, 0.352 and y « 0.034. An enlargement of this long orbit is shown in fig. 8.) Is 
it the KAM surface around this island chain that will determine the asymptotic behavior of 
/t? In this picture, the asymptotic behavior is determined by the islands at a finite depth in 
the islands-around-islands hierarchy. This may be false. Perhaps as t is increased, we must 
look at deeper and deeper levels of the hierarchy. Such is the view taken by Chirikov in ref. 
13. We may also have to jump between branches of the hierarchy as t increases. (This is 
supported by the observation that the three longest orbits seen for K = 0.1 are each stuck 
close to different island chains.) 

The absence of any idea as to the asymptotic behavior of Cr makes it virtually impos- 
sible either to obtain numerically an upper bound on the diffusion coefficient or to establish 
that it is unbounded. Such questions can probably only be answered when this problem is 
better imderstood analytically. In any case, it may be necessary to take a critical look at 
extraneous effects which act as an extrinsic noise source and which will cut off the very long 
correlations. 

Finally, what is a measure of these islands? This gives the overall importance of this 
phenomenon. We already have Sinai's estimate as to the number of islands that will be 
created. What is still needed is some estimate of their size and of the interval in parameter 
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spacc in which they exist . Numerieal evidence suggests that these both decrease rapidly with 
the period. The latter quantity would enable us to say what is the measure of parameter 
space on which we expect to see islands. 
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Appendix A: Properties of the periodic mapping. 

In order to obtain correct trapping statistics, we would like to uniformly sample the 
trajectories the area outside the islands of Q and measure how long each trajectory spends 
close to the islands. We achieve this goal by looking at long trajectories of Q*. In this case, 
"close" is determined by Xmin and Xmax- However, fig. 2 shows that the area outside the 
islands is not sampled uniformly. There are holes in it corresponding to the spurious islands. 
Suppose that the order of a particular set of islands is n. Particles in this set of islands 
will have looped around the torus at least once in the x or y direction after n iterations of 
Q*. (These islands are accelerator modes.) A trajectory of length n within these islands 
is made up of one or more segments. It cannot be part of a longer segment. By omitting 
the trajectories within the islands we are therefore undercounting the segments with lengths 
less than or equal to n. 

We can express this in a more quantitative way. The total area of phase space in Q* 
is bo = L^. We define bi to be that part of 60 which is not occupied by the islands of Q 
and &2 to be the area of the stochastic component of Q*, i.e., that part of bo which is not 
occupied by the islands of Q* . The difference between these two, 61 — 62 gives the area of 
the spurious islands. To include the effect of the spurious islands, we define /(* by 

ft T-ft + ~^ 9t, (Al) 

where gt is the contribution from the regions of phase space occupied by the spurious 
islands. If gt is normalized like ft to have a first moment of unity, then /* has the same 
normalization. The factors multiplying ft and gt are just the relative areas of those regions 
of phase to which these trapping statistics apply. Now gt specifies the lengths of orbit 
segments in the spurious islands, which must be less than or equal to Timax the maximum 
period length for these islands. Thus gt is zero for t > rimax- If we wish to calculate Ct, we 
need only know the trapping statistics for t > t. So as long as r > rimax, we may ignore gt 
and use 

Ol 

These spurious islands have one other interesting effect. They introduce correlations 
between the lengths of successive segments. Consider an orbit close to such an island chain 
of period 3. Then, in the simplest case, we expect to see a succession of segments of length 
3. Such correlations are obviously a peculiarity of Q*. In the general mapping, where the 
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particlc spends some time in the stoeliastic sea before returning to the neighborhood covered 
by Q* , the trapped segments will be sampled randomly. 

We next examine the effect of the definition of "close," i.e., at the effect of changing Xmin 
and Xmax- Consider the effect of increasing the box size by decreasing ccmin and increasing 
aimax- This has two effects: new segments which never entered the original box appear; 
and all the original segments may be extended both forwards and backwards. As long as 
the original box was large enough, the new segments will all be short (since they never get 
close to the islands). By the same argument, the extension of the orbit segments will be 
small. Because of the extremely rapid departure of orbits away from the islands, we can, 
for instance double, the edge of the box and only increase the average length of the orbits 
segments by less than one. So the main effect of increasing the box size will be to increase 
the number of short segments. As in the case of the accelerator modes, we can quantify 
this. If we change Xmin to a;^;^ and a;inax to x'^^y., then, for large enough t, the new trapping 
statistics fl are given by 

b'^n = b2fu (A3) 
where h'2 is the area of the stochastic component in the new box. 

Lastly, we address the problem of estimating the trapping statistics from an orbit of 
length T. Consider an infinitely long time record consisting of segments of various lengths 
t — 1,2^ . . .. Pick a random record of length T and call Nt the expected number of segments 
of length t. (We do not count partial segments at the ends.) Then 



for i > T, 

(T + 1 - t)ft for < t < T. 



This may be proved by induction. When T = 0, we have A^t = as required. When we 
extend the record from T — 1 to T, the only way a new segment of length t can be included 
is if the particle is at the end of such a segment at time T. Because the probability of this 
happening is /j, Ni increases by fi for t <T. 

Thus we divide actual number of orbit segments by {T + 1 — t) to obtain an unbiased 
estimate of the trapping statistics, ft- The proof given above does not depend on the 
segments appearing in a random order. This is important because, as we have seen, there 
will be correlations between successive segments in Q*. One problem with the method we 
use for measuring Nt is that the beginning of the record is not arbitrary, because initially 
we have {X, Y) uniformly distributed on the line X = 0. If T is large, this is not expected 
to affect the results very much. 



Appendix B: Numerical methods used. 

In order to make the reduction of x and y in Q* to their base intervals as cheap 
as possible, the L x L square is reduced to a unit square with a zero origin. Thus new 
coordinates X and Y are defined by 

x = Xmin + LX, y=-\L + LY. 

In these coordinates Q* becomes 

N: Yt^ Yt-i + G{Xt-i) (mod 1), Xt = Xt-i + Yt-\ (mod 1), 

where 

G{X) = aX'' -hX + c, a = 2L, b = -Ax^i^, c=2{xl,^- K)/L. 
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Wc will always choose .T,„i„ < —\fK, so that the stable fixed point at {x,y) = (— \/^, 0) is 
included in the base intervals. Thus the quantities a, 6, and c are all positive. 

N is implemented using fixed-point arithmetic on the Cray-1. This should be compared 
with the use of integer mappings by Rannou [14]. In both cases "exact" (in a sense to 
be defined later) numerical mappings can be defined. The principal difference is in the 
coarseness of the grid on which the mapping is defined. The finest grid that Rannou used 
was 800 X 800. In contrast, our grid is about 2**® x 2"'^, so that its resolution is close to 
that obtained with floating point numbers. The use of fixed-point arithmetic also results in 
faster performance because the extraction of the integer and fractional parts of a number 
just correspond to masking operations. 

First we define the notation for fixed-point numbers. We view the full 64-bit word as 
representing a number in twos-complement binary notation with the "binary" point placed 
before the last 48 bits. Addition and subtraction are performed by the integer add and 
subtract instructions. Multiplication of fixed-point numbers is accomplished by the floating- 
point multiply instruction provided the two numbers lie in the range [0, 1). Multiplication 
of a number by 2™ where m is a non-negative integer is accomplished by shifting the word 
to the left m places (with zeros introduced on the right). 

The function G is then computed by writing it as 

G{X) = 2"^{AX^ -BX + C), (A, B, C) = 2-™(a, 6, c) 

where m is the smallest non-negative integer such that A, B, and C all lie in [0,1). X 
also lies in this range so that all the multiplications in G may be carried out by computing 
AX^ — BX + C with the add, subtract, and multiply instructions and shifting the result 
m places to the left. Since the low m bits of G are always zero, then, assuming that the 
low m bits of both X and Y are initially zero, they will remain zero. Thus there will be 
n = (^2'^s-m^2 pQii^^g accessible in the unit square. (Typically m = 2 and n k, 10^^.) 

Wc are now ready to give the procedure for gathering the trapping statistics ft . 

(1) Read in K, Xmin, a;max, X^., T, q. (The trapping statistics will be collected for the qth 

power of the mapping. T is the total number of iterations of N~'i. X^ defines the size 
of the randomizing zone.) 

(2) Initialize 5^ ^ for < i < 150, t ^ 0, r ^ T, {X,Y) ^ (0,0). Initialize the 
random number generator with a "random" seed. {Bi are the bins used for collected 
the trapping statistics. The counter t gives the length of the orbit segment so far. The 
number of steps still to do is given by r.) 

(3) Set L ^ <— 2L, b ■* 4a;min, c ^ 2(0;^;^ — K)/L. Calculate m. A, B, C, 

according to the the prescription given above. 

(4) Set T ^ T - 1. If r < stop. 

(5) Initialize P ^ false. (P is set to true when the trajectory leaves the base square.) Set 

t^t + 1. 

(6) Set {X, Y) ^ R{X, Y) where the mapping R is given below. 

(7) Set {X,Y) ^ N^'^{X,Y) where the mapping A'^"^ is given by iterating the mapping 
N^^ given below q times. 

(8) If P is false go to step 4. 

(9) Record new orbit segment by setting i <— b{t), Bi Bi + 1. Re-initialize f •*— 0. Go to 

step 4. 

The mapping R is defined by 
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(6.1) liX> Xr return. 

(6.2) Set P ^ true, Y Y + U where f/ is a uniformly distributed random variable 
(mod 1). Return. 

The mapping A''"^ is defined by 

(7.1) Set X ^ X - r + i. If [XJ ^ 0, set P ^ true. Set X ^ X - \_X\ . ([XJ is the largest 
integer satisfying X > [X\ .) 

(7.2) Set y ^ y - G(X) where G(X) is calculated as outlined above. If [FJ ^ 0, set 
P ^ true. Set y ^ y - [yj . Return. 

The bin number b{t) is defined by 

{t for t < 50, 

[SOlogioiJ for50<^<10^ 
150 for 10^ < t. 

Finally the trapping statistics ft are obtained by 

'Bi/{T+l-i) fori<50, 

Fi= { B,/(r+l-10('+i/2)/30) 



together with 



|"10(»+i)/30l _ |"ioV30-| 
ft = Fi/q\ 



otherwise, 



where i = b{[t/q\). 

The mapping R randomizes the Y positions of the particles in a thin zone (0 < X < X,.) 
on the left edge of the square. {Xr is chosen to be small, typically 0.02.) As long as the 
randomizing zone is sufficiently far removed from the islands, R has no effect on the length 
of orbits which are trapped for moderately long times. R serves two purposes. When it 
is first applied, it picks a random initial condition on the line X = 0. It also stops the 
trajectory from being periodic either because of the spurious accelerator modes or because 
a stochastic trajectory has eventually returned to its starting position. This results in a 
more rapid sampling of phase space and so better statistics are obtained. 

Since the accumulation of the statistics Bi is quite expensive compared with iterating 
the map, we look instead at the qth power of the mapping. We count the orbit segment as 
having ended if the orbit wrapped around the torus in any of the q iterations. This leads 
to some uncertainty in ft for f ~ g. However, we are most interested in very long trapping 
times for which t ^ q. 

Full advantage is taken of the vectorization capabilities of the Cray-1 by following M (a 
multiple of 64) orbits and averaging the results. Finally, the critical parts of the algorithm 
were handed-coded in Cal, the Cray Assembly Language. The execution time is approx- 
imately (360 + 72>q)T ns per particle. This is about 6 times faster than a straightforward 
implementation in FORTRAN using floating-point arithmetic. The most time-consuming 
computations were those of fig. 6 where M = 1600 orbits of length qT = 2 x 10^ were used, 
a total of 3.2 x 10^^ iterations of N. This consumed about 65 hours of CPU time on the 
Cray-1. 

Aside from the computation of G{X), all the calculations are exact. N is therefore area- 
preserving, or, since only a discrete set of numbers may be represented on the computer, it 
is better to say that it is one-to-one. A trajectory may be run backwards by 

TV"* = J-^TV* J 
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where J is the involution ( is the identity) 

J: Xt=Xt-i-Yt-i + \ (modi), Ft = -^t-i (modi) 

which may be calculated exactly using fixed-point arithmetic. 

Because of the discreteness of the number system on a computer, all the trajectories 
of a numerical mapping are eventually periodic. Because iV is a one-to-one mapping, the 
whole trajectory will be periodic. There will be no initial transient. Presumably the typical 
period length of an orbit in the stochastic part of phase space is very long. Rannou [14] gives 
^(n-l- 1) as the average length of a trajectory in a random permutation of n points. (This is 
to be compared with the period length for an orbit in a random function on n points which 
is 7rn/8 + ^ approximately [15]. The distinction between a function and a permutation 
is that a permutation is a one-to-one mapping of the n points onto themselves, while a 
function is a many-to-one mapping of the points into themselves. It is not clear whether the 
stochastic orbits of floating-point mappings are closer to those of random permutations or 
random functions.) However, N is not a random mapping because it possesses symmetries. 
In the case of the standard mapping this reduces the expected length of stochastic orbits to 
0(\/n) [14]. Although this number is still large (about 10"'^'^), we expect there to be a large 
number of trajectories with shorter periods. Examples of such trajectories are those in the 
accelerator modes which are restricted to a small fraction of phase space. By composing the 
mapping with i?, we give the trajectory a random step whenever it returns to its starting 
point (and usually before that). We therefore ensure that a particle stays in these short 
periodic orbits for at most q periods. (Of course the random number generator is an example 
of a mapping. But it is specially chosen to have a very long period of 2^'^, so the random 
number generator repeats itself after about lO^^g iterations of N .) 



Appendix C: Correlation function for Markov model. 

In section 6, we presented a model for the motion of a stochastic trajectory when an 
island is present. We wish to make that model more precise so that we can accurately assess 
the influence of the island on the correlation function. We do this by modeling phase space 
with a finite number of points and writing down a Markov transition matrix P for evolution 
of the distribution function. 

Suppose there are a total of N points. In the simplest case, there is just one trapped 
segment of length m. The stochastic sea consists of the other n = N — m points, one of 
which is special in that it is the pre-image of one of the trapped points. For example if 
N = 7, m = 3, n = 4, then P is given by 




















1\ 


1 























1 























1 

4 


1 

4 


1 

4 


1 

4 











1 

4 


1 

4 


1 

4 


1 

4 











1 

4 


1 

4 


1 

4 


1 

4 










1 

4 


1 

4 


1 

4 


1 
4 





(CI) 



The entry in the jth row and ith column is the probability that a particle at point i at time 
t is at point j at time t + 1. Thus we see that points 1-3 constitute the trapped segment, 
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whilc points 4-8 are the stochastic points with point 8 being the pre-image of one of the 
trapped points. 

In the more general case, there may be M trapped segments, whose lengths are 
for 1 < fc < M. Suppose the members of the fcth trapped segment are Tfc = {p^j^\p^^\ • • • , 
pj™'"''}, and that the pre-image of p^."^-* is qk. Wc will also define = {qk\ U T^, T = (Jj, T^, 
T' — UfcT^, S = T'^, S' = T"=, where the superscript c denotes the complement, and 
m = ruk- Thus T is the set of trapped points, T' is T with the pre-image points added, 
S is the set of stochastic points, S' is the S with the pre-image points removed, and m is the 
total number of trapped points. It is easy to generalize (CI) to incorporate the additional 
trapped segments. 

If V{t) is a column vector giving the distribution function at time t, then 

V{t) = P*F(0). 

As t oo, P* approaches a matrix C each of whose entries is the constant 1/iV. It is 
convenient to define Q such that Q = P — C. It is easy to show that, for f > 0, Q* = P* — C 
so that limj^oo Q* = 0. Since so many of the entries in Q* are the same, the computational 
effort involved in doing the matrix multiplication is the same as for multiplying (m + 1) x 

(m -I- 1) matrices. 

The correlation function for a function h on phase space is given by 

C(t) = h • P^ • h/N, (C2) 

where h is a vector giving the value of h for each point. This parallels the definition (13) 
given in section 6. As in that section, we will assume that J2 = 0- ^^^^ case, we have 
h • C • h = 0, so that we can replace P by Q in this definition. 

Let us consider first the case where there is a single trapped segment, M = 1. Suppose 
that hj is 1 for j e T, for j = q\, and takes on arbitrary values (consistent with the 
requirement that Yl^j = 0) for j G S". In the example given in (CI), wc might choose 
h = [1,1,1,-1,1, —3, 0]. (We pick h to be zero for the pre-image point, so that there is no 
bias for this point.) In the notation of section 6, we have Ai = N, Bi = m, f* = Stm/'m, 
Ci,(T) = (m - t)/N for r < m, Vi, = \m^/N. 

We have numerically determined C(r) for this case using (C2) with m = 10, 20, and 
50 and N varying between 10m and lO^m. In particular, we considered e(r) = (C(r) — 
Cis(T))/Cis(0). We found that for r > 0, |e(r)| < m/N = Bi/Ai. After an initial transient, 
e(r) settles down to a decaying periodic function, which oscillates about zero, whose half- 
period is m, and which decays as n""^/"* approximately. Thus for N ^ m, the corrections 
to Cis(r) are very small. 

It is, perhaps, a little surprising that this model predicts anti-correlations where C(r) < 
0. This is most pronounced for m < t < 2m. It is, however, a reflection of a real effect 
in mappings. Consider an orbit in the stochastic sea of a general mapping. If there is an 
island present, the trajectory can be stuck close to the island for very long times. However, 
because of ergodicity, the average time the orbit spends close to the island is equal to the 
relative areas of the stochastic components Bi/Ai. Thus, to compensate for the stickiness, 
an orbit far from the island must have some difficulty coming close to the island. This 
leads to the anti-correlations observed. (These anti-correlations may be rather difficult to 
measure directly because the presence of trapped segments of so many different lengths will 
average out the oscillations in C.) A fairly accurate picture of the situation is that the 
stochastic region close to the island is surrounded by a box with a very small opening in 
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it. A particle' inside the box has to bounce around a lot before finding the opening and 
escaping. Similarly, once the particle is outside the box (in the main part of the stochastic 
sea), it cannot easily find the opening in order to come close to the island once again. 

Because Cis(T) = for r > m, whereas c{t) is finite for all times, we still need to verify 
that the integrated effect of e(T) (which gives the corrections to Vis) is small. We define the 
diffusion coefficient in the same way as (6). This gives 

oo 

P = h • (i| + ^ Q") • h/7V = h • (i| + Q • (I - Q)-i) • h/N. 

T=l 

We are able to write the sum in this way because Q'" converges to zero. The the matrix op- 
erations were performed for several examples with 1 or 2 trapped segments using Macsyma 
[16]. The results show that in general V may be written as 

^ = ^(E^^+EEE^a)- (C3) 

This result is exact and only requires that hj = 0. The first term in parentheses is just 
the contribution from C(0) for the stochastic points, while the second term gives island con- 
tribution. In the one-segment case discussed above, the island contribution clearly reduces 
to ^w? /N which exactly equals Ag. Thus the oscillations in e(r) are such that its sum is 
precisely zero. 

For the more general case, where there is an arbitrary number of trapped segments 
of various lengths, there are two requirements: hj should be a constant ft-o for j T and 

nikhq^. should be zero. (This last condition says that the pre-image points should be 
sufficiently evenly spread over phase space.) From section 6, we have /j* = '^i-Stmk/'f^j 
Vis = \hQ/N^f,m\, while from (C3) we obtain 

^ = ^(E^'+^oE-^'^)- m 

^j£S k ^ 

Again, the second term exactly gives Ag. Alternatively, we could declare that the pre-image 
points are part of the trapped segments. (This conflicts with the picture given in section 
6 because it would allow one trapped segment to follow immediately after another with no 
intervening free segment.) Letting hj = ho for j G T' and defining = + 1, (C4) would 
be modified by replacing S by S' and irik by mj^,. 
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Fig. 2. (a) Stochastic trajectories for periodic quadratic mapping Q* for K = 0.1. (b) An 
enlargement of a portion of (a). Here Xmin + VK = —0.4, x^ax + VK = 0.6, Xr = 0.02. 
(a) was produced by plotting every lOOOth point of 64 orbits each of length 10^ and (b) by 
plotting every 100th point of 64 orbits of length 2 x 10*. 
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Fig. 3. The trapping statistics ft for (a) K = -1 0"^, (b) K = 0, and (c) K = 10"^. 

In each case x^i^ + ^ max(i^', 0) = —0.5, x,„ax + \/ niax(/r, 0) = 0.5, = 0.02. For (a) 
M = 128, q = 1, qT = 5 X 10^; for (b) and (c) M = 512, q = 2, qT = 10^. In (b) and (c) 
the straight hnes give the t~'^ and exp(— O.lf) behaviors predicted in section 4. 




Fig. 4. Trajectories, for the Hamiltonian with (a) K = (7) and with (b) K small and 
positive (11). In both figures H takes on equally spaced values with an increment of |. 
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Fig. 6. (a) The trapping statistics ft for K = 0.1. (b), (c), and (d) show Pt, Cr, and 
rflogCr/dlogr for the same case. Here Xmin + = —0.4, x^ax + = 0.6, X^. = 0.02. 
The data for ft was obtained in 3 pieces; for 1 < t < 10^, M = 128, g = 1, qT = 5 x 10^; 
for 10^ < t < 10^, M = 256, g = 10, gT = 5 x 10^; for 10^ < t, M = 1600, q = 1000, 
qT = 2x 10^. 
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Fig. 8. An enlargement of a long orbit segment in fig. 6(a). The bottom left part of the 
picture is occupied by a member of a chain of 138th-order islands. 



